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1.   INTRODUCTION 

This  report  outlines  the  design  of  a  general  purpose  graphic 
interactive  modeling  and  simulation  system  in  which  the  user  can  define  the 
primitives  of  the  system  in  the  first  stage,  and  then  use  them  to  model  his 
problem  in  the  second.   This  is  similar  to  the  compiler-compiler  concept  in 
that  a  language  is  used  to  define  a  further  language.   The  first  language 
allows  the  definition  of  the  elements  while  the  second  consists  of  the  use  of 
these  elements .   It  is  unlike  the  compiler-compiler  in  the  requirement  that 
the  first  language  must  also  be  heavily  user  oriented,  so  that  every  user  will 
be  readily  able  to  define  his  own  elements  during  the  development  of  his  model, 
although  he  may  make  use  of  libraries  of  fairly  common  elements. 

The  use  of  graphics  during  all  stages  of  system  use  is  stressed  as 
an  important  part  of  the  interactive  quality  of  the  system,  although  it  is 
certainly  true  that  with  suitable  input  and  output  conventions,  the  system 
could  be  operated  from  a  typewriter  terminal  or  through  batch. 

In  any  modeling  system  the  user  is  provided  with  a  number  of 
modeling  blocks  (elements)  which  follow  certain  laws  (that  is,  obey  certain 
equations)  and  he  may  connect  these  in  various  ways.   These  connections  may 
be  of  a  topological  nature,  that  is,  elements  may  be  joined  together  at 
certain  connection  points  (called  terminals  below)  or  may  satisfy  some  global 
relation  (for  example,  have  a  common  ambient  temperature).   Existing  modeling 
systems  provide  predefined  elements,  although  -they  may  allow  the  user  to 

define  additional  elements  of  a  similar  kind.   Examples  of  such  systems  exist 
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in  simulation  (e.g.  the  IBM  CSMP  system   )  in  which  the  elements  represent 

typical  analog  computer  elements,  and  in  network  analysis  (e.g.  the  Electronic 

Circuit  Analysis  Program  ECAP).   These  systems  are  tailored  for  particular 

classes  of  problems,  and  cannot  be  used  on  larger  classes  without  a  great  deal 


of  analysis  on  the  part  of  the  user.   The  purpose  of  the  proposed 
system  is  to  make  it  possible  for  the  user  to  define  any  type  of 
element  and  the  way  in  which  it  can  interact  with  other  elements. 
It  should  not  "be  restricted  to  single  variable  elements  as  are 
simulator  elements,  or  to  two  terminal  elements  as  are  electrical 
network  elements.   Thus,  the  elements  might  represent  atoms 
reacting  dynamically  with  a  number  of  near  neighbors  in  a  lattice, 
cars  passing  through  an  intersection,  or  complex  electrical  network 
elements  with  many  terminals . 

A  typical  simulation  program  will  allow  the  user  to 
join  a  number  of  elements  together  in  a  network  such  as 


In  this  diagram,  element  1  is  an  inverter,  elements  2  and  3  are 
integrators,  while  element  k   is  a  gain.   These  elements  are 
characterized  by  the  fact  that  the  input  and  output  "wires"  each 
"carry"  the  value  of  a  single  real  number,  and  that  they  are 
directional  elements;  that  is,  they  have  distinct  inputs  and  out- 
puts.  Thus,  the  input  to  element  1  carries  the  value  y,  so  its 
output  carries  the  value  -y.   Since  element  2  is  an  integrator, 
its  output  z  is  given  by  the  expression 

t 


■  pi  +/0<-y>« 


The  P  is  a  parameter  that  can  be  set  as  an  initial  value.   The  output 

of  element  h   is  P  times  its  input,  or  -P~y.   The  output  of  element  3  is 

the  integral  of  the  sum  of  its  inputs  with  the  initial  value  of  P9  so  we 

have 

.t 
0 
On  differentiating  twice  we  get 


y  =  P2  +J  (z  -P3y)dt 


ix.  =  _y  .  P  dv_ 

dt2    y    3  dt 

y(o)  =  p2 

^L(O)  =  P  -  P  P 
dtv  ;    1    2r3 

Thus,  the  particular  network  can  be  used  for  integrating  a  linear  homo- 
geneous second  order  differential  equation.   In  addition  to  these  elements, 
many  others  which  combine  the  inputs  in  both  continuous  and  discontinuous 
fashions  are  available  and  most  systems  will  allow  the  user  to  define  some 
new  elements  by  means  of  available  procedure  oriented  language  subroutines. 
In  all  cases,  however,  the  elements  are  uni-directional;  that  is,  they  have 
inputs  and  a  single  output  each  of  which  carries  only  a  single  value.   The 
output  of  any  element  may  be  connected  to  any  number  of  input  elements,  but 
not  to  another  output  element,  since  then  the  number  carried  on  the  "wire" 
would  be  doubly  defined. 

These  requirements  also  lead  to  other  types  of  restrictions,  the 
most  important  being  that  no  loop  may  exist  that  does  not  contain  an  element 
that  causes  a  delay.   A  loop  is  any  closed  path  such  as  occurs  in  Figure  1 
going  from  element  1  to  2  to  3  and  back  to  1.   A  delay  can  be  caused  by  a 


delay  element  that  states  that  the  output  is  the  value  of  the  input 
at  a  given  earlier  time  or  an  integrator.  This  type  of  restriction 
insures  that  the  situation  shown  below  cannot  occur. 


p  p 

*-*© -d> 


This  network  defines  y  to  be  P  P  y,  which  means  that  y  is 
zero  unless  P-,Pp  =  1.   Most  digital  simulators  retain  this  restric- 
tion because  it  simplifies  the  method  of  solution — on  an  analog 
computer  such  a  connection  is  not  practical  because  of  inherent 
delays  in  the  physical  realization  of  the  element. 

Systems  have  been  produced  to  interface  graphic  terminals 
with  simulator  programs.   See,  for  example,  Baskin    .   The  user  is 
presented  with  a  set  of  predefined  elements,  and  is  allowed  to  define 
additional  elements  that  are  of  the  same  general  class. 

Network  analysis  programs  use  two  terminal  elements  that 
are  bi-directional;  that  is,  there  is  no  distinction  between  input 
and  output.   Whereas  the  nodes  of  a  network  for  a  simulation  system 
were  governed  by  the  fairly  simple  requirement  that  a  single  output 
signal  must  be  connected,  and  that  this  determines  the  value  of  all 
connected  inputs,  the  nodes  of  the  electrical  network  are  governed 
by  Kirchoff's  laws — namely,  that  the  voltage  of  all  lines  connected  to 
the  node  is  the  same  and  that  the  sum  of  the  currents  entering  the 
node  from  all  connected  elements  is  zero.   The  former  can  be  replaced 
by  the  statement  that  the  sum  of  the  voltage  drops  across  elements 


in  any  closed  loop  is  zero.  These  facts  are  used  to  take  the  equations 
that  govern  each  element  and  derive  a  system  of  differential  equations. 
A  knowledge  of  the  element  equations  is  built  into  the  network  analysis 
systems  and  cannot  be  changed. 

In  a  similar  manner,  structural  networks  consisting  of  members 
(columns  and  beams,  etc.)  are  connected  together  at  nodes.   The 
individual  members  are  governed  by  equations  (Hookes  law  relating  the 
elastic  "stretch"  to  the  applied  force,  strain  and  stress  to  give  them 
their  correct  terms,  for  example)  while  the  ends  of  the  beam  are  charac- 
terized by  twelve  physical  quantities,  the  six  components  of  deflection 
and  the  six  components  of  force t    The  nodes  imply  the  equation  that  the 
sum  of  all  of  the  forces  acting  on  a  node  is  zero,  component  by  component, 
and  that  the  deflection  of  all  members  connected  to  a  node  is  the  same. 
Thus,  there  is  a  direct  analogy  between  current  and  force  and  between 
voltage  and  deflection.   In  fact,  the  same  methods  of  analysis  can  be  used 
for  the  two  problems . 

In  the  analog  computer,  a  restriction  that  no  loops  could  be 
specified  that  did  not  include  a  delay  was  enforced.   No  such  restriction 
is  reasonable  in  the  case  of  electrical  or  structural  networks;  consequently, 
the  network  analysis  programs  contain  techniques  to  analyze  such  loops  and  to 
calculate  their  overall  effect.   The  network  analysis  programs  can  perform 
this  extra  service  because  they  "know"  the  type  of  behavior  possible  in  each 
element,  and  can,  therefore,  apply  known  techniques  for  solving  the  resulting 
equations.   The  simulation  program  users  are  allowed  to  specify  arbitrary 
elements  by  means  of  computer  programs  so  it  is  difficult  to  determine  when 
a  closed  loop  will  have  a  solution. 


The  proposed  system  will  allow  arbitrary  elements  to 
be  specified  and  used  in  networks,  so  it  will  contain  standard 
simulators  and  network  analyzers  as  a  subset.   Naturally,  the 
generality  will  lead  to  a  loss  of  speed  and  it  is  appropriate  to 
ask  at  this  point  how  this  will  be  overcome  and  how  the  system 
will  be  able  to  handle  the  more  complex  networks  that  are  not 
allowed  in  existing  systems.   Because  the  system  is  not  yet 
fully  implemented,  this  can  only  be  partically  answered.   New 
numerical  techniques  for  the  solution  of  ordinary  differential 
equations    will  be  used.   These  have  produced  speed  increases 
of  several  orders  of  magnitude  for  one  class  of  problems 
encountered  in  simulation.   The  use  of  sparse  matrix  techniques 
and  an  initial  symbolic  processing  of  the  equations  can  replace 
repeated  numerical  passes  over  parts  of  the  equations  by  a  single 
symbolic  pass.   This  is  discussed  in  some  detail  in  Section  k. 
The  matrix  techniques  will  also  be  needed  for  the  detection  of 
incorrectly  specified  networks.   Because  the  usual  restrictions 
of  simulation  and  network  analysis  systems  are  dropped,  it  is 
possible  for  the  user  to  give  networks  with  many  solutions  or  none. 
It  is  possible  to  detect  this  in  almost  all  cases  with  techniques 
given  in  Section  h.      The  unsolved  problem  is  how  to  guarantee  that 
all  cases  are  detected,  and  more  important,  how  to  relate  detected 
under  or  over  specifications  to  the  original  user  input.   When  a 
user  of  an  analog  computer  simulator  fails  to  connect  an  output  to 
one  of  the  nodes,  the  simulator  can  tell  him  of  the  error  in  those 
words.   When  the  system  is  represented  by  a  large  set  of  differential 
equations,  the  system  will  only  detect  that  there  are  more  variables 


than  equations.   The  solution  will  require  some  restrictions  on  element 
specification,  but  it  is  not  clear  what  they  are  at  this  time. 

Section  2  will  discuss  the  graphical  interface  between  the 
user  and  the  system,  Section  3  the  numerical  algorithms  used  and 
Section  5  indicates  what  the  current  state  of  the  art  is  by  describing 
the  capabilities  of  the  system  that  can  (and  are)  being  implemented, 
without  further  theoretical  advances. 
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2.   GRAPHICAL  SPECIFICATION 

An  interactive  graphic  display  will  be  used  to  specify  the 
elements  used  and  the  network  connections.   The  use  of  graphics  during 
the  basic  element  definition  stage  is  for  mnemonic  value  only,  whereas 
the  graphical  description  of  the  element  connections  provides  a  more 
concise  description  of  the  network  than  most  other  techniques.   Graphics 
will  also  have  an  obvious  application  during  the  data  input  and  output' 
stages . 

In  very  complex  networks  in  which  some  sections  behave  almost 
independently  of  the  others,  it  may  be  helpful  to  allow  the  user  to 
point  out  subareas  that  can  be  separately  analyzed. 


2a.   DEFINING  BASIC  ELEMENTS 

An  element  will  be  represented  by  a  mnemonic  picture  which  is 
drawn  by  the  user.   The  element  contains  a  number  of  parts: 

1.   An  arbitrary  line  drawing. 

This  is  the  mnemonic  for  the  element  which  is  used  for 
identification  purposes  by  the  user.  It  could  consist 
of  the  drawing 


for  example.   There  is  no  significance  to  any  of  the  lines, 

so  the  picture  is  retained  within  the  graphic  section  entirely. 

2.  Arbitrary  character  strings. 

These   are  part   of  the   line    drawing  and  have   only  mnemonic 
value.      The   figure   above  might   also   contain  the   string 
"l  WATT",  presumably  to  identify   it   as    a  one  watt   component. 
This   naming  should  not  be   confused  with  parameter  specifica- 
tions  described  below. 

3.  Labeled  terminal  points. 

Any  of  the  points  in  the  line  drawing  may  be  labeled.   We  will 
use  labels  1,  2,...  etc.   Each  label  must  be  unique.   The 
drawing  above  could  have  its  left  and  right  ends  labeled  1  and 
2,  respectively.   These  terminals  will  be  the  connection  points 
of  the  element . 
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h.      Descriptive  equations. 

These  specify  the  action  of  the  element.   They  can  involve 
variables  defined  on  the  terminals  and  variables  internal 
to  the  element.   If  the  figure  above  refers  to  a  resistor 
which  depends  on  temperature,  the  equations  might  be 

E(2)  =  E(l)  +  (RO  +  R1*T)*I(2) 

Kl)  =  -K2) 

T'    =  (E(2)  -  E(1))*I(2)*K  -  C*(T  -  TO) 

Here  the  resistance  is  assumed  to  be  a  linear  function 
of  the  element  temperature  T,  which  itself  is  determined 
by  a  first  order  differential  equation  involving  the 
voltage  drop  (E(2)  -  E(l))  and  the  current  1(2).   TO  is 
the  ambient  temperature,  while  C,  K,  RO,  and  Rl  are 
constants  for  the  device.   Prime  denotes  differentiation. 

5.   Parameters. 

In  most  applications,  elements  such  as  the  one  specified 
will  be  used  repeatedly  with  different  values  of  the 
device  constants.   Therefore,  these  must  be  parameters 
and  be  given  values  each  time  they  are  used.   Consequently, 
at  the  specification  stage  the  user  can  also  place  labels 
on  the  figure  corresponding  to  the  symbols  used  in  the 
equations.   At  the  time  of  use,  these  will  be  replaced  by 
numeric  values  or  expressions  involving  other  variables. 
Thus,  the  complete  specification  for  this  example  is 
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"1  WATT" 


©-AAAJ  Vr© 


C  K  RO  Rl 

E(2)    =   E(l)    +    (RO  +   R1*T)*I(2) 

KD    =   -K2) 

Tf         =    (E(2)    -   E(1))*I(2)*K   -   C*(T   -   TO) 

Variables  were  assumed  to  "be  associated  with  each  of  the  nodes. 
In  the  examble  above,  voltage  E  and  current  in  I  was  assumed  on  each  of 
the  two  nodes.   In  a  general  problem,  there  can  be  one  or  many  variables 
associated  with  each  node.   A  loose  jointed  three  dimensional  structure 
has  six,  for  example,  while  an  analog  computer  has  one.   Thus,  the  system 
must  be  told  how  many  variables  appear  on  each  terminal  and  what  their 
names  are.   In  a  later  step,  terminals  of  several  elements  will  be  joined 
together  in  a  network  node.   Rules  are  needed  for  generating  relations 
between  the  variables  on  the  connected  terminals.   In  all  problems  studied 
by  this  writer,  only  two  types  of  rules  have  been  necessary:   Kirchoff's 
current  and  voltage  laws.   Thus,  we  will  identify  variables  as  being  I  type 
or  E  type  depending  on  which  rule  they  follow. 

It  may  also  be  possible  for  different  types  of  terminals  that  may 

not  be  interconnected  to  occur.   (This  could  arise  because  of  different 

signal  levels,  for  example.)   In  that  case  the  user  would  like  the  option 

of  specifying  which  types  of  terminals  are  to  be  used.   Therefore,  before  any 

elements  are  specified,  the  user  will  have  to  answer  the  following  questions: 

HOW  MANY  NODE  TYPES  WILL  YOU  USE? 

NAME  ALL  E  TYPE  VARIABLES  ON  NODE  TYPE  1. 

NAME  ALL  I  TYPE  VARIABLES  ON  NODE  TYPE  1. 
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(The  latter  two  questions  are  repeated  for  each  node  type.) 
In  the  example  above  the  answers  would  be  1,  E,  and  I. 

The  system  will  then  proceed  through  the  sequence  of 
questions i 

DO  YOU  WISH  TO  SPECIFY  A  NEW  ELEMENT? 
If  the  answer  is  yes,  the  opportunity  to  make  an  arbitrary  line  drawing 
is  provided.   (it  is  also  possible  to  include  character  strings.)   An 
escape  from  this  will  lead  to  the  question: 

DO  YOU  WISH  TO  INDICATE  A  TERMINAL  POINT? 
A  yes  answer  will  yield  the  opportunity  to  indicate  a  point  in  the 
figure.   It  is  labeled  by  the  next  available  integer,  starting  at  1, 
and  the  user  is  asked 

WHAT  IS  ITS  NODE  TYPE? 
(unless  only  one  node. type  was  specified). 

When  all  terminal  points  have  been  specified,  the  user 
is  asked: 

DO  YOU  WISH  TO  SPECIFY  PARAMETERS? 
If  yes ,  he  indicates  their  positions  by  light-pen  and  inputs  their 
names.   Finally  he  is  asked  to 

PLEASE  SPECIFY  THE  ELEMENT  EQUATIONS 
He  will  type  them  in  the  manner  indicated  in  the  above  example. 
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2b.   SPECIFYING  THE  NETWORK 

When  all  desired  elements  have  been  specified,  the  user  will 
be  in  a  position  to  use  them  in  a  network.   Each  element  will  appear  as 
a  light-button  (menu  item)  on  the  display.   Instances  of  the  elements 
may  be  selected  and  positioned  on  the  display.   Terminals  of  elements 
may  then  be  connected  to  form  nodes  of  the  network.   If  more  than  one 
node  type  is  specified,  the  system  will  restrict  interconnections  to 
terminals  of  the  same  type.   As  instances  are  connected,  the  system  will 
be  producing  systems  of  equations  which  describe  the  network. 

As  an  instance  of  an  element  is  selected,  the  user  is  asked  to 
specify  the  values  of  each  parameter.   This  can  be  by  numeric  input 
(assigning  the  value  of  a  resistor,  for  example),  or  the  parameter  could 
be  given  another  symbolic  value  for  later  specification.   (An  element 
could  depend  on  its  temperature  which  is  later  defined  to  be  the  temper- 
ature of  another  element,  for  example.) 
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2c.   RESTRICTIONS  ON  THE  ELEMENTS  AND  THE  NETWORK 

The  result  of  a  set  of  element  descriptions  and  a  network 
which  uses  them  will  be  a  set  of  algebraic  and  differential  equations 
It  is  evident  that  it  must  be  possible  to  solve  these  if  the  network 
is  to  have  any  meaning.   (There  are  cases  in  which  the  solution  is 
not  unique  but  for  which  such  non-uniqueness  is  unimportant.   These 
will  be  discussed  in  Section  3. ) 

Consequently,  it  must  be  possible  to  determine  whether  or 
not  the  system  can  be  solved.   This  is  probably  the  most  severe 
problem  introduced  by  the  generality  of  the  processor.   There  is 
naturally  a  trade-off  between  the  difficulty  of  determining  solv- 
ability and  the  freedom  allowed  the  user.   In  existing  systems,  the 
number  of  variables  allowed  on  terminals,  the  type  of  relations 
which  can  be  given  and  connection  restrictions  may  be  built  in.   A 
simulator  system,  for  example,  allows  one  variable  of  E  type  on 
each  terminal.   It  may  have  many  terminals  and  some  terminals  may  be 
distinguished  as  output  terminals.   One  and  only  one  output  terminal 
must  be  connected  to  each  node  of  the  network.   Furthermore,  many 
simulation  systems  restrict  the  network  connections  to  avoid 
"unnatural"  feedback  connections — those  loops  that  contain  no  delay 
elements.   Thus  if  the  circle  represents  a  gain  element  whose  output 
is  g  times  the  input  defined  by 


Q^-^Q © 


E(2)  =  g*E(l) 
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the  network 


is  not  one  that  would  normally  be  allowed  on  an  analog  computer,  although 
algebraically  it  serves  to  define  E(l)  =  E(2)  =  0.   (if,  however,  the 
gain  were  unity,  it  would  not  define  the  value  of  E  so  in  an  algebraic 
sense  it  would  retain  its  initial  value.   It  would  be  unwise  to  use  such 
a  subnetwork  in  an  analog  computer  program! )   On  the  otherhand,  electrical 
networks  restrict  the  user  to  two  terminal  elements  (usually)  with  two 
variables  on  each  terminal.   It  is  also  known  that  the  I  type  variable 
sums  to  zero  for  the  element  and  that  the  E  type  variable  enters  only  via 
the  difference  between  the  values  on  the  two  terminals.   There  are  no 
distinguished  output  terminals,  although  there  are  restrictions  against 
connecting  voltage  generators  in  parallel  or  current  generators  in  series. 

(Parenthetically,  it  should  be  noted  that  simulator  elements 
can  be  viewed  as  a  set  of  two  terminal  elements.   Thus,  the  element 


C  =  a*A  +  3*B 


can  be  viewed  as   three   two-terminal  elements 
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where  the  terminals  1,  2,  and  3  are  "grounded"  in  use  and  the 
equations  are 

C  _  E(3)  =  a*(A  -  E(l))  +  $*(B  -  E(2)) 
The  I  variable  for  each  terminal  can  be  ignored.  ) 

If  no  restrictions  are  placed  on  the  definition  of  elements, 
it  is  necessary  to  wait  until  the  complete  network  has  been  specified, 
collect  all  relevant  equations,  and  then  attempt  to  determine 
whether  or  not  they  are  singular.   By  placing  some  restrictions  on 
the  definition  of  new  elements,  the  process  of  determining  non- 
singularity  can  be  reduced.   The  system  requires  that  the  restrictions 
given  below  are  adhered  to  in  the  definition  of  elements  and  their 
connection  in  networks .   Alone  they  are  not  sufficient  to  ensure  that 
the  equations  are  non-singular,  but  they  ease  the  analysis  of  the 
network. 

Restrictions  on  Elements 


forms 


or 


The   equations    for  the   internal  variables   u  may  take  the 


(i)      u'    =   f(all  variables) 


(ii)    u     =   a(all  other  variables) 
where   the  prime    denotes    differentiation  with   respect  to  time.      There 
must  be   exactly  one   such  equation   for  each  internal  variable.      The 
terminal  variables  v  must   satisfy  equations   of  the   form 

v  =   a(all   other  variables) 
where   each  v  may  only  appear  on  the   lefthand  side   of  one   such  equation 
(but  may  not   appear  on   any). 


IT 


An  option  in  which  the  equality  sign  is  replaced  by  the 
left  arrow  («-)  assignment  operation  allows  the  user  to  indicate  an 
element  which  "forces"  the  value  of  a  terminal  variable.   Only  one 
such  terminal  may  be  connected  to  any  node.   Thus,  the  equality  sign 
is  used  in  the  strict  sense  rather  than  the  procedural  language 
assignment  sense. 

Network  Restrictions 


Restrictions  which  prevent  multiple  assignment  of  variables 
on  a  node  by  the  left  arrow  are  obvious.  Similarly,  equations  of  the 
form 

E(l)  =  0 
are  equivalent  to  assignments,  but  can  easily  be  detected  and  so 
treated.   More  difficult  are  equations  of  the  form 

E(l)  =  E(2)  +  1 
(that  is,  voltage  generators).   Two  of  these  may  not  be  connected  in 
parallel,  but  recognition  of  this  and  more  complicated  forms  have  to 
be  left  to  the  analysis  stage  described  in  Section  h. 
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2d.   COMPOUND  ELEMENT  DEFINITION 

It  is  desirable  to  be  able  to  use  previously  defined 
elements  in  the  definition  of  new  elements.  The  steps  in  this 
process  are  as  follows : 

1.  Construct  a  network  from  existing  elements* 

2.  Identify  nodes  in  this  network  as  terminals  of 
of  the  compound  element. 

3.  Specify  any  additional  equations  connected  with 
the  element.   These  may  include  new  element 
variables  which  could  be  parameters  in  some  of 
the  constituent  elements.   For  example,  TO,  the 
ambient  temperature  in  the  resistor  example  above, 
may  vary  for  the  component  element  due  to  heating 
effects  of  all  of  the  elements.   It  could  also  be 
specified  by  differential  equations. 

k.  Replace  the  drawing  of  the  network  by  a  new  line 
drawing  that  serves  to  identify  the  element,  and 
indicate  the  terminal  points  on  the  line  drawing. 

5.   Specify  parameters  of  the  compound  element,  and 
indicate  where  they  are  to  be  placed. 

The  new  drawing  will  then  be  available  as  an  element  for 
the  construction  of  further  networks. 


1') 


2e.   DATA  INPUT 

Prior  to  the  simulation  of  the  network,  values  must  be  provided 
for  any  varialbes  as  yet  unspecified  and  certain  nodes  may  have  to  be  set 
to  given  input  values.   Nodes  can  always  be  set  by  defining  "voltage 
source"  elements  with  a  given  time  dependent  behavior,  so  the  problem  is 
equivalent  to  one  of  specifying  arbitrary  function  values  for  some 
variables .   Thus ,  we  may  have  the  element 


E(l)  =  E(2)  +  F 

where  F  is    a  variable   and  we  wish  to   assign  F   a  value.      If  it   is   a  simple 
function,    it   can  be   expressed  algebraically,    for  example 

"F  =    3*SIN(T)" 

However,  we  wish  to   allow  more   general  inputs.      It  must  be 
possible   to   specify   an   arbitrary   function  F(Y)   by  means   of  a  table   or  a 
graph,    and  to  use   this    function  in  the   definition  of  a  variable. 
Consequently,  we   allow  the   input   of  tables   or  single  valued  graphical 
functions   and  use   the   information  to   determine   an   approximation  to  the 
function.      It   is   not   obvious  what   approximation  should  be   used,   but   a 
Spline   function  would  have   the   advantage   of  derivative   continuity.      It  may 
also  be   desirable  to   allow  discontinuous    functions   although  these   slow  the 
simulation  process   down  because   of  numerical  problems  with  high  order  methods 

When   all  parameters  have  been  specified,    it  is   necessary  to 
determine   the   initial   conditions.      The   user  may  wish  to   specify  some,    others 
can  then  be   determined  by  performing  a  simulation  with   constant   inputs   until 
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the  system  settles  (if  it  is  stable).   An  integration  with  time  as  the 
independent  variable  can  then  be  used  to  follow  the  system  behavior. 
Values  of  terminal  or  element  variables  or  expressions  containing  them 
can  be  specified  as  output  which  can  be  obtained  in  tabular  or  graphic 
form. 

It  is  obviously  desirable  to  be  able  to  repeat  analysis  for 
different  parameter  and  initial  values,  so  the  final  step  is  to' allow 
statements  in  a  conventional  programming  language  to  be  used  to  write 
loops  and  to  call  the  simulation  system. 
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3.   EFFICIENCY,  NUMERICAL  CONSIDERATIONS 

When  a  complete  network  has  been  specified,  the  equations 
representing  each  of  the  elements  can  be  collected.   These  equations 
will  involve  the  variables  on  the  elements  and  on  their  terminals.   The 
equations  consist  of  a  set  of  linear  and  nonlinear  algebraic  and 
differential  equations.   They  will  be  written  in  the  form 

ui     =  fi(uj'TJ'"rt)  (3-1) 

Ti     =  «L(YVVt)  (3-2) 

R. .w.  =  (P. .u.  +  Q. .v.)  (3.3) 

where  the  prime  represents  differentiation  with  respect  to  time,  the 
U  are  the  set  of  variables  defined  on  the  elements  which  satisfies 

differential  equations,  the  v.  are  the  set  of  variables  that  appear 

either  on  elements  or  on  terminals  that  are  specified  by  nonlinear 
equations  and  the  w.  are  the  remaining  variables.   Summation  is 

presumed  to  occur  on  the  j  subscript  in  the  third  set  of  equations. 
If  the  network  is  improperly  specified,  then  these  linear 
equations  may  have  no  solution  or  may  possess  a  non unique  solution. 
There  may  be  more  or  fewer  equations  than  there  are  unknowns.   (in 
other  words,  R  may  or  may  not  be  square  and  even  if  square  may  not  have 
an  inverse.)   However,  even  if  the  equations  are  properly  specified,  there 
may  not  be  a  unique  solution.   If,  for  example,  the  variables  on  the  nodes 
are  actually  current  and  voltage  and  all  elements  only  specify  the  voltage 
drop  across  the  element  (there  is  no  element  that  specifies  the  actual 
voltage  of  one  terminal),  then  the  value  of  the  voltage  of  at  least  one 
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node  will  be  undetermined.   All  other  nodes  may  be  given  as  a  level 
relative  to  that.   Similarly,  if  all  current  sums  for  elements  are 
zero,  at  least  one  of  the  Kirchoff  current  laws  for  the  nodes  is 
redundant.   Thus,  although  the  matrix  R  may  be  square,  it  could  be 
singular  and  the  equations  have  many  solutions.   For  this  section 
we  will  assume  that  R  is  square  and  nonsingular.   The  difficulties 
will  be  discussed  in  Section  k. 
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3a.   STIFF  DIFFERENTIAL  EQUATIONS  AND  NONLINEAR  ALGEBRAIC  EQUATIONS 

The  resulting  system  of  differential  equations  is  frequently- 
stiff.   This  means  that  there  are  widely  separated  time  constants  present. 
This  causes  most  common  methods  to  be  extremely  slow,  using  steps  which 
are  of  the  order  of  the  shortest  time  constant  present  even  though  the 
components  due  to  these  time  constants  have  been  damped  out  for  a  long 
time.   Recently  some  new  methods  have  been  used  very  successfully  on 
such  problems  .     They  are  a  class  of  methods  of  varying  orders  of  the 
multistep  type,  but  can  be  organized  so  that  both  the  step  and  order  can 
be  chosen  automatically  in  order  to  make  the  step  as  large  as  possible 
while  controlling  the  single  step  error.   They  require  no  special  starting 
procedure  because  it  is  part  of  the  step  and  order  control  mechanism. 
These  methods  will  be  used  in  the  system. 

If  a  system  of  the  form  (3.1)  is  given — we  will  ignore  the  other 
variables  for  the  moment — it  is  necessary  to  calculate  an  approximation  to 

the  Jacobian  J  =  9f./8u.  in  order  to  solve  the  corrector  formula.   The 

i   J 

iteration  is  given  by 


(m+1)  _  _.(m)  ,  ,_   v  rx  ,  ^  Tl-1    r..(m)   _   ,..(m) 
Xik 


■•*£  +  \*  i»m    Ki-w,  <«;:•  'i 


J 


where  u   is  the  calculated  approximation  to 


,k-l  k-1 
h    d   u. 

l 


(k-l)!dtk_1 


and  h  is  the  time  step.   These  approximations  are  used  to  estimate 
(predict)  the  value  of  u.,  at  the  next  time  step.   The  coefficients  b, 

IK  i 
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determine   the   method.      The    corrector  iteration   actually   amounts  to   a 
Newton-Raphson  iteration  if  the  Jacobian  J  is   re-evaluated  at  each 
iteration,    otherwise   it   is    a  chord  method.      As   long  as   it   converges, 
it   converges   to   a  solution  of  the   corrector,    so  it   is   not  necessary 
to  have   a  precise   value   of  J.      The   system  of  nonlinear  algebraic 
equations    also  has   to  be   solved  at  the   same   time.      These    can   also  be 
handled  by   a  Newton  or   a  chord  method  in   a  similar  way.      Suppose   that 
a  good  first   approximation  to  v        (defined  in  the   same  way  as   u      )    is 

known.      Then  the   iteration 

J  J    Ij 

converges  to  the  set  of  numbers  v..  in  which  v.,  will  be  the  solution 

lk  ll 

of  the   algebraic  equations,   while   the   remaining  v.,    will  be   approx- 
imations to  the   scaled  derivatives.      (They  will  be   in  error  by  the 
first  neglected  derivative   after  sufficient   steps.)      Again  the 
Jacobian  need  not  be   re-evaluated  at  every  step   of  the   iteration, 
because   it   only   affects   the   rate   of  convergence,   not   the   final  solution. 
The   reason   for   calculating  the   additional   derivatives   is   that  very 
accurate   first   approximations   to  the  value   at   the  next  time   step   can  be 
obtained  so  by  prediction  that  very  few  iterations    are  needed.      (One 
or  two   is  typical.)      The    constants    c      determine  the   order  of  the  method, 

K. 

that  is,  how  many  of  the  scaled  derivatives  are  needed.  A  combination 
method  in  which  the  corrector  equation  and  the  algebraic  equations  are 
solved  at  the  same  time  is  given  by  calculating  the  matrix 
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1   +   b     h3f/3u 


Z  = 


cl  9£./3iL 


b  +  h  3f/3v 


I  +   c,    9g_/3y_ 


-1 


and  forming  the   correction   factors 


d  =   Z 


4m)    -hf(uU),   v(m),   t) 


-1 


/    (m)         (m)      ,  v 


where  u  is  the  vector  formed  from  u._,  while  vn  =  v  is  the  vector 
— id  id  — l   — 

formed  from  v.,.   Products  of  the  elements  of  d_  with  the  vectors  b, 

and  c   can  then  be  added  to  the  matrices  u.,   and  v..   to  form  the 
k  lk       ik 

(m+l)-st  iterate. 

It  was  ass-umed  above  that  the  vectors  f  and  g_  are  functions 

of  u  and  v  only.   They  are  also  functions  of  w  and  this  must  be  taken 

into  account  in  calculating  Z.   Since  we  assumed  that  R  is  nonsingular 

we  can  write 


w  =  R-1(Pu  +  Qv) 


(3.U) 


We  will  write 
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^=  A,  3^=B,  3^=  C 


3u 


3v 


3w 


^  =  D,  *_£=   E,  ^=  F 
3u      3v      3w 


We  may  then  write 


I  +  b  h(A  +  CR  P) 


Z  = 


cx(D  +  FR  1P) 


b  h(B  +  CR  -"-Q) 


-1 


I  +  c.,(E  +  FR  \) 


We  also  need  to  calculate  the  functions  f  and  g_  for  a  given  u,  v, 
and  t.   To  do  this  we  must  evaluate  (3.^0.   The  linearity  of  u  and 
v  in  this  is  not  important,  so  we  could  have  allowed  (3-3)  to  be 
written 


Rw  =  ^(u,  v) 


(3.5) 


where 


P  =  ^  and  Q  =  ^ 


3u 


3v 
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3b.   SPARSE  MATRICES 

Typical  networks  contain  many  elements  and  even  more  equations 
involving  only  a  few  variables.   Consequently,  the  work  involved  in 
calculating  Z  and  evaluating  (3.^)  may  be  the  limiting  factor  in  speed. 
For  this  reason  we  use  sparse  matrix  techniques.   Since  the  multiplica- 
tion of  Z  by  a  vector  occurs  in  every  iteration  step,  the  objective  is 
to  compile  code  to  perform  this  as  rapidly  as  possible,  and  never 
actually  generate  the  elements  of  Z  explicitly.   The  value  of  the  Jacobian 
may  only  be  recalculated  periodically  (when  the  step  size  changes,  for 
example)  so  that  the  calculation  that  is  dependent  only  on  the  Jacobian 
should  be  done  separately  and  not  repeated.   This  is  the  subject  of  the 
next  section. 
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3c.   GENERATING  THE  LINEAR  EQUATIONS 

The  linear  equations  can  be  solved  directly  by  matrix 
inversion  (using  a  sparse  matrix  technique)  so  it  is  desirable  to 
place  as  many  equations  in  the  form  (3.3)  as  possible.   Since  we 
allowed  them  to  take  the  more  general  form  (3.5)  we  must  give  an 
algorithm  for  determining  which  set  of  variables  is  w.   The  set  u 
can  be  determined  immediately  by  examining  the  element  equations.- 
The  remaining  element  equations  and  node  equations  must  first  be 
scanned  to  remove  the  trivial  relations,  and  then  scanned  to  generate 
the  nonlinear  equations.   When  a  nonlinear  equation  of  the  form 

x  =  x(all  variables) 
is  given  we  must  check  to  see  if  all  the  variables  that  appear 
nonlinearly  in  the  equation  have  already  appeared  on  the  lefthand 
side  of  the  differential  or  nonlinear  equations  already  generated. 
If  not,  one  of  these  variables  can  be  used,  say  v.,  and  the  equation 

J 

written  as 


v.   =  v.    +  x(all  variables)   -  x  (3.6) 

J    J 


(In  fact,  it  is  not  necessary  to  add  the  v.  to  both  sides,  but  is  is 

J 

helpful  for  expository  purposes  to  make  it  similar  in  appearance  to 
the  differential  equations.)   If  all  of  the  variables  appearing  non- 
linearly have  already  appeared  on  the  lefthand  side,  the  equation  can 
be  saved  as  a  candidate  for  a  "linear"  equation.   When  the  scan  has 
been  completed,  a  set  of  such  candidates  remains.   (Parenthetically, 

we  note  that  the  choice  of  which  variable  v.  to  place  on  the  left  of 

J 

(3.6)  may  affect  the  number  of  equations  left  as  possible  linear 
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equations.   Thus,  if  we  have  the  system 


0  =  x  +  x  *x  +  1 


0  =  x  +  x  **2  +  x 


0  =  xi  +  l/x2  -  x3 


and  we   choose   to  write   the   first   as 


X2  =   x2  +  Xl  +  x2*x3  +   -1 


the  remaining  two  can  be  written  in  the  form 


Xl  +  X3  =  X2**2 


whereas  had  we  started  with 


x~  =  x  +  x  +  X  *X_  -  1 


we  would  also  have  called  the  second  nonlinear  and  written 


x2  =  x2  +  x1  +  x2**2  +  x3 


leaving  one  linear  equation.   The  solution  technique  will  give  the  same 
answers,  but  the  speed  may  be  different.   One  of  the  questions  yet  to  be 
investigated  is  that  of  optimizing  the  assignment.) 

The  remaining  equations  can  be  written  in  the  form 

Rw  =  q(u,  v) 

but  R  may  be  singular.   It  is  therefore  necessary  to  remove  some  of  the 
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singularities .   Consider  the  part  of  a  network  below 

[1 


representing  two  inductors  in  a  series .   These  will  eventually 
reduce  to  the  equations 

1(1)'  =  (E(l)  -  E(2)/L1 

1(2)'  =  (E(2)  -  E(3))/L2 
and 

1(1)  -  1(2)  =  0 

These  serve  to  determine  (E(l)  -  E(3))and  to  determine  E(2)  in 
terms  of  E(l)  and  E(3)  (independent  of  l(l)  and  1(2)).   Unfortun- 
ately the  scheme  given  above  would  place  l(l)  and  1(2)  in  the  set 
u  and  E(2)  in  the  set  w.   The  third  equation  would  be  called  linear 
so  its  row  of  R  would  be  zero,  making  R  singular.   If  the  third 
equation  is  placed  in  the  "nonlinear"  class  and  E(2)  placed  in  the 
class  v  by  writing 

E(2)  =  E(2)  +  1(1)  -  1(2) 
the  numerical  technique  will  work  for  nonzero  step  sizes  h. 

Consequently,  rows  and  columns  of  R  must  be  removed.   This 
will  have  to  be  detected  during  the  symbolic  inversion  process 
discussed  in  the  next  section. 
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An  alternative  technique  that  can  be  used  in  this  case  is  to 
recognize  the  problem  during  the  scan  to  remove  trivial  equations. 
This  would  try  and  replace  1(2)  by  l(l)  and  would  lead  to  two  differential 
equations  for  l(l).   The  second  could  be  converted  to  the  algebraic 
equation 

(E(l)  -  E(2))/L1  =  (E(2)  -  E(3))/L2 
which  will  give  one  differential  equation  and  one  linear  equation  for 
this  section. 

This  latter  technique  cannot  eliminate  all  of  the  singularities 
in  R.   For  example,  an  element  between  L  and  L  above  could  have  given 

a  nonlinear  relation  between  l(l)  and  1(2).   Thus,  we  must  be  prepared 
to  detect  singularities  in  R  during  the  analysis. 
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k.      NETWORK  ANALYSIS—SYMBOLIC  MANIPULATION 

A  number  of  steps  have  to  be  taken  symbolically  because  a 
key  to  the  system  is  the  ability  to  perform  as  much  of  the  matrix 
inversion  prior  to  the  integration  as  possible  and  to  be  able  to 
detect  singularities  in  the  system.   Once  numbers  have  been  inserted 
we  cannot  guarantee  that  zeros  arise  from  singularities  rather  than 
cancellation, 

The  first  step  is  to  develop  the  Jacobian  matrix  for  all 

equations.   This  requires  a  partial  differentiation  as  performed  in 

[h] 
ODESSY    .   It  is  then  generally  simple  to  detect  constants  which 

indicate  linear  terms.   (it  may  not  be  possible  to  find  all  such 

constants  because  full  automatic  simplification  is  not  possible.) 

This  information  is  used  to  generate  the  tentative  set  of 

linear  equations,  and  hence  the  matrix  R.   The  next  step  is  to  try 

and  invert  R  "symbolically."   If  singularities  are  found,  they  must 

be  removed  by  taking  out  the  appropriate  rows  and  columns  and 

placing  them  in  the  set  of  nonlinear  equations.   There  is  an  except 

tion  to  this  :   If  during  the  reduction  of  R  a  zero  is  encountered, 

and  the  righthand  side  arising  from  p(u,  v)  is  also  identically  zero, 

it  means  that  one  of  the  variables  in  w  can  be  arbitararily 

specified  and  that  the  zero  equation  can  be  ignored  completely.   This 

occurs,  for  example,  when  all  elements  involve  only  a  voltage 

difference  and  all  current  sums  are  zero.   One  or  more  of  the  Kirchoff 

current  equations  can  be  eliminated,  and  one  or  more  of  the  voltages 

can  be  chosen  arbitrarily.   In  this  case,  the  equation  need  not  be 

placed  in  the  nonlinear  group. 
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Ha.   SYMBOLIC  SPARSE  MATRIX  INVERSION 

We  are  given  a  matrix  R  whose  elements  may  be  zero,  small 
integers,  constants,  or  symbolic  expressions.   The  problem  is  to 
solve  the  system 

RX  =  B 

where  R  is   K  x  N,   X  and  B  is   N   x  M,   M  _>  1.      Typically,   the  number  of 
nonzero  elements   of  R  is   small    (2N  to   5N),   N  is   large    (>  100),    and 
M  small    (often  l).      The   elements   of  B  may  also  take   one   of  the   four 
forms    cited.      The  values   of  the   "variables"   appearing  in  the  symbolic 
expressions  will  be   specified,    at  which  time  the   value   of  X  is   needed. 
The   goal  is   to   do   as   much   "preprocessing"   as   possible.      Two   situations 
need  to  be   considered  in  the   analysis.      In  the   first   the  values   of  B 
will  change  very  frequently,  while  the  values   of  R  will  only   change 
infrequently.      In  the   second,  we  need  to   calcualte   a  third  M  x  M 
matrix  of  the   form 

D  =   CR_1B 

where  C  is  M  x  N,  and  in  this  case,  R  may  be  singular  (so  that  R  B  can 
be  interpreted  as  any  solution  of  RX  =  B,  B  being  assumed  to  have  the 
appropriate  properties  for  such  a  solution  to  exist).   If  rows  of  R  lead 
to  singularities  such  that  no  solution  exists ,  these  rows  must  be 
identified  and  removed  from  R.   Since  D  itself  will  later  take  part  in  a 
similar  problem,  the  solution  for  D  should  be  as  non-numeric  as  possible— 
by  this  is  meant  that  all  elements  that  will  be  of  the  form  of  small 
integers  or  simple  symbolic  expressions  should  be  carried  that  way.   This 
will  be  clearer  after  reading  the  example  following. 


3k 


EXAMPLE 


Suppose  we  have 

R(I,1)   =   1     -1     0  0     S+T 

R(I,2)    =0        L      0  0      -L 

R      -R      1  0      0 

1        0      0  10 

R(I,5)   =0        10  0      1 


and 


B(l)  =  1 
B(2)        =   0 

C 

1 
B(5)        =   D 

We  will  identify  R(6,j)  with  B(j)  so  that  R(l,J)  is  a  5  x  6  matrix. 
The  Gaus  elimination  back  substitution  method  (without  concern  for 
best  pivots)  leads  to  the  following  steps  to  convert  the  elements  of 
B  to  the  elements  of  the  solution  X: 

-R*R(I,1)  +  R(I,3).+  R(I,3)    I  =  1,...,  6 

Note  that  the  symbolic  operation  for  1=2  should  give 

(-R)*(-l)  +  (-R) 
=  R*l  -  R 
=  R  -  R 
=  0 
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-R(l,l)  +  R(I,U)  +  R(I,U) 

Note  that  this  makes  B(U)  =  -1  +  1  =  0  since  we  can  do 
exact  arithmetic  on  small  integers. 

These  two  steps  reduce  the  5  x  6  R  to 

1-10      0  S+T  1 

0      L        0      0  -L  0 

0      0        10  -R*(S+T)  C-R 

0      10      1  -(S+T)  0 

0      10      0  1  D 

Divide  R(I,2)  by  L 

Note  that  this  makes  R(5,2)    into   -1  if  symbolic   arithmetic 
is  performed.      Note   that  B(2)    is    zero,    so   is   not   changed. 

-R(I,2)   +  R(I,U)   ->  R(l,U) 

Note   that  B(2)    is    zero,    so  B(U)    is   not   changed. 
-R(I,2)   +  R(I,5)   ■+  R(I,5) 

Note   that  B(5)    is  not   changed  since  B(2)    is    zero,    and  that 
R(5,5)   becomes   -(-l)   +1=2,   exactly. 

Divide  R(l,5)   by  2 
These   steps   reduce   R  to 

1-10     0     S+T  1 

0        10     0-1  0 

0        0      10      -R*(S+T)  C-R 

0        0      0      1      -(S+T)+l  0 

0        0      0      0      1  D/2 
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The  back  substitution  process   is    as   follows 

(S+T-1)*R(I,5)  +  R(I,^)  ■+  R(I,U) 
R*(S+T)*R(I,5)  +  R(I,3)  ■>  R(I,3) 
R(I,5)  +  R(I,2)   -*■  R(I,2) 

-(S+T)*R(I,5)  +  R(I,1)  -  R(I,1) 
R(I,2)  +  R(I,1)   -*  R(I,1) 

which  reduces   R  to 

10      0      0      0  l-(S+T)*D/2   +D/2 

0      10      0      0  D/2 

0      0      10      0  C-R  +  R*(S+T)*D/2 

0      0      0      10  (S+T-l)*D/2 

0      0      0      0      1  D/2 

The  solution  X  appears  in  the  last  column. 

If  this  were  to  be  done  for  large  systems ,  it  is  evident  that  the 
symbolic  forms  would  get  unreasonably  complex,  even  if  optimum 
simplification  were  to  be  performed.   Since  the  objective  is  to  get 
fast  code  which  means  that  common  expressions  should  be  factored 
out,  it  is  sensible  to  remove  them  at  an  early  stage  in  the  process. 
Because,  in  this  example,  R(l,5)  =  S+T  is  already  complex,  it  might 
be  better  (it  is  in  this  case)  to  rename  it  R.. ,_  which  is  equivalent 

to  saying  that,  "We  have  a  way  of  calculating  this  element  but  we  do 
not  know  anything  about  its  relation  to  other  elements  in  the  matrix." 
This  may  mean  that  if  S+T  appeared  elsewhere  we  could  lose  the 
possibility  of  cancellation,  but  if  expressions  become  very  complex 
the  problem  of  automatic  symbolic  simplification  is  not  tractable  in 
a  reasonable  amount  of  computer  time.   Note  that  saying  that,  "We 
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have   a  way  of  calculating  the  number,"   is   implied  if  we   say,    "We 
have   compiled  code   to  produce   its   value   in   a  known   location   in 
memory. " 

We   do  not    consider  the   fact  that  we   are,    in  some   cases, 
going  to  want  to   form  R"1  many  times  where   the   values   of  L,   R,    S, 
and  T  are   fixed   (so  that   R"1  is    fixed)   but  where   C   and  D   (and  hence 
B)    change.      This    can  be   done  by  the  two   sets   of  steps  below.      The 
first   and  second  gets   must  be   done   each  time   that  L,   R,   S,    and  T 
change;    only  the   second  set  need  be   done   if  C  and  D   change.      These 
steps    are   obtained  by  going  back  over  those   steps   given   above   and 
writing  down  only  those   that   result   in   "complex"   expressions;   that 
is,   ones  we   do  not  wish  to   retain  symbolically.      In  the   steps   below, 
all  names   refer  to  the   contents   of  locations   in  memory  which   do  not 
need  to  be   stored  in   any  particular  order. 

SET  1  SET  2 


1. 

S+T 

-> 

R!5 

2. 

"R*B15 

-> 

R35 

C-R  +  B3 

3. 

"w 

■+■ 

% 

5. 
6. 

7. 

8. 

9. 


D/2  -*■  B 

-VR35*B5  +  B3 

1-R15*B5  ♦  Bl 
B1"B2  +  \ 
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Thus,   we   use    a  multiplication   and  two   additions    (or 
subtractions)   to   form  the  new  "operator  inverse   of  A,"  that   is   a 
sequence   of  steps   S  which  when   applied  to  B  give  R     B.      (i.e. 
S(B)   =  R     B).      The   sequence   S   is   given  in  SET  2,    and  takes   one 
division,   three  multiplications. 

Two  problems  have  been  ignored.      They  relate  to  the 
selection  of  the  pivot  element.      In  the  example   above   it  was 
always   taken  to  be   the  next    diagonal  element.      The   choice   of 
pivot   affects  both  the   retention  of   zero  elements    (which  is 
obviously   desirable   in  order  to   reduce  the   amount   of  work)    and 
the   numerical   accuracy — in  the  worst   case   the   diagonal  element 
could  be    zero   so   that   the   division  of  the   row  by  the   diagonal 
element   is  not  possible.      To   see   the  effect   of  bad  selection  of 
pivots   on  the    zeros,    consider  the   5x5  matrix 

1  x  x  x  x 

x  1  0  0  0 

x  0  1  0  0 

x  0  0  1  0 

x  0  0  0  1 

where  x  represents  any  nonzero  matrix. 

If  the  (l,l)  element  is  used  as  pivot  the  first  time,  the 
resulting  matrix  will  be  of  the  form 

1  x  x  x  x 

0  x  x  x  x 

0  x  x  x  x 

0  x  x  x  x 

0  x  x  x  x 
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resulting  in  the  loss  of  all  of  the  initial  zeros.   (This  is  the 
worst  case!).   If,  instead,  the  pivots  had  been  taken  in  the  order 
(2,2),  (3,3),  (h,h),    (5,5),  and  (l,l),  no  zeros  would  have  been  lost 
in  the  process . 

Numerical  error  requirements  dictate  that  the  pivot  element 
may  not  be  "too  small."   The  exact  meaning  of  too  small  is  complex 
and  beyond  the  scope  of  this  report,  but  in  most  programs  the  largest 
element  in  the  chosen  column  is  used.   (Since  any  row,  representing 
an  equation,  can  be  scaled  by  any  nonzero  quantity,  it  is  obvious  that 
each  row  must  be  normalized  in  some  manner  to  apply  this  criterion.) 
Our  desire  is  to  perform  most  of  the  inversion  symbolically  before  the 
actual  values  are  known  so  it  is  not  possible  to  make  such  a  choice 
in  general. 

Consequently,  we  must  make  some  reasonable  assumptions  during 
the  symbolic  steps,  and  possibly  note  some  tests  that  should  be  made 
during  the  subsequent  numerical  steps  to  verify  the  reasonableness  of 
the  operations.   We  started  with  either  numerical  values  (integers  or 
constants)  or  symbolic  expressions.   Since  we  cannot  tell  much  about  the 
value  of  expressions,  we  must  assume  that  they  could  take  the  value  zero 
(although  it  is  evident  that  expressions  such  as  X2  +  1  are  bounded 
away  from  zero),  and  that  we  must  try  and  use  elements  with  numerical 
values  as  pivots  first.   Pivots  with  the  value  1  or  -1  are  most  desirable 
because  no  division  is  necessary.   (The  distinction  between  small  integers 
and  constants  is  simply  the  question  of  certainly  of  value.   If  we 
perform  the  symbolic  manipulation  (-2)*L  +  2*L  where  -2  and  2  are  "small 
integers,"  we  know  we  get  zero,  whereas  if  we  perform  (-2.0)*L  +  2.0*L 
where  -2.0  and  2.0  are  constants  that  have  been  calculated  and  happen  to 


Uo 

be  the  same,  we  only  know  that  each  represents  a  number  whose 
magnitude  is  in  the  range  2.0  +  e,  consequently,  the  expression 
could  be  in  the  range  +  2eL.   It  may  prove  necessary  to  use  this 
fact  in  automatic  analysis . ) 

The  two  goals  of  pivot  selection  may  lead  to  conflicting 
criteria:   The  largest  element  in  a  column  may  cause  the  greatest 
destruction  of  zero  elements.   Therefore,  in  choosing  the  pivot,  we 
will  first  concentrate  on  avoiding  zero  loss,  then  on  numerical 
accuracy.   (This  may  appear  backwards  since  the  end  result  has  first 
to  be  accurate  then  to  be  rapidly  obtained,  but  there  is  some 
flexibility  in  pivot  choice  in  many  problems — only  a  very  bad  choice 
will  cause  unacceptable  errors — whereas  bad  pivot  selection  often 
increases  the  number  of  nonzeros ,  and  hence  operations,  drastically.) 

The  choice  of  pivot  elements  has  received  study  for  sparse 
matrices  (ones  with  many  zeros)  by  a  number  of  people.   It  does  not 
appear  to  be  possible  to  find  the  best  strategy  without  examining  all 
permutations  of  columns  (which  is  not  practical  for  large  N),  but  a 

number  of  heuristics  have  been  tried  experimentally. 

[3] 
Dantzig,  et.al.    tested  ten  different  rules  on  linear 

programming  problems.   Although  no  one  rule  is  uniformally  better,  one 

was  recommended  for  their  class  of  problems.   We  will  use  it  until  tests 

have  been  made  to  determine  better  rules . 

The  rule  is  first  to  examine  each  column  to  find  out  how  many 

nonzero  elements  it  contains.   The  subset  { C}  of  those  columns  with  the 

minimum  number  of  nonzero  elements  is  chosen.   Of  all  rows  with  nonzero 

elements  in  the  columns  (c),  the  first  row  A  with  the  minimum  number  of 

nonzero  elements  is  chosen.   The  pivot  element  is  the  first  common 

nonzero  element  of  A  and  {C}. 
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We  will  modify  the  last  sentence  to  "the  pivot  element  is  the 
most  simple  element  common  to  A  and  {C},  where  most  simple  means  an 
order  in  which  1  and  -1  are  simplest,  followed  by  small  integers,  variables, 
constants  reasonably  different  from  zero,  and  lastly  symbolic  expressions." 
Note  that  after  the  first  pivot  has  been  chosen  and  the  necessary  row 
operations  performed,  the  next  choice  should  be  made  on  the  basis  of  the 
of  reduced  matrix,  ignoring  the  rows  already  chosen  as  pivot  rows. 

We  also  wish  to  bias  the  selection  in  favor  of  rows  with  zeros 
in  the  B  matrix.   Consequently,  in  choosing  the  row  when  there  are 
several  rows  with  the  same  number  of  elements,  we  favor  the  first  with  a 
zero  row  in  B. 

We  will  illustrate  the  mechanization  of  this  technique  on  our 
earlier  example.   To  simpify  the  operations  of  searching  for  the  best 
rows  and  columns  we  will  initially  count  the  number  of  nonzero  elements 
in  rows  and  columns : 


NON   ZERO 

COUNT 

3 

k 

1 

1 

3 

1 

3 

1 

-1 

0 

0 

X 

1 

2 

0 

L 

0 

0 

-L 

0 

3 

R 

-R 

1 

0 

0 

c 

2 

1 

0 

0 

1 

0 

1 

2 

0 

1 

0 

0 

1 

D 

x  stands  for  an  arbitrary  number. 

The  first  step  gives  columns  3  and  k   with  one  element,  Row  h 
is  chosen,  so  the  first  pivot  is  the  (U,U)  element.   No  operations  are 
required  on  the  matrix,  but  not  the  element  in  the  (l+,l)  position  is  no 
longer  considered  so  the  form  is 
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STEP 
NUMBER 


SC 


NONZERO 
COUNT  ZC 


2   1+   1    13 


SR     ZR 


3 
2 
3 

2 

2 


1 

-1 

0 

0 

X 

1 

0 

L 

0 

0 

-L 

0 

R 

-R 

1 

0 

0 

c 

1 

0 

0 

1 

0 

1 

0 

1 

0 

0 

1 

D 

The  order  in  which  rows  and  columns  are  chosen  is  also 
specified  in  vectors  SR  and  SC.   The  next  pivot  is  (3,3)  which  also 
requires  no  row  operations,  but  ZC(l)  is  reduced  to  1,  ZC(2)  to  3, 
while  SR(3)  =  SC(3)  =  2.   Consequently,  column  1  is  the  next 
choice,  and  there  is  only  one  nonzero  element  in  that  column  in 
the  rows  not  yet  handled,  so  it  is  the  pivot;  that  is,  the  (l,l) 
element.   Now  the  matrix  is 


SR     ZR 


3      3 

2 

2      3 

1      2 

2 


SC 

3 

- 

2 

1 

- 

ZC 

1 

2 

1 

1 

1 

-1 

0 

0 

X 

1 

0 

L 

0 

0 

-L 

0 

R 

-R 

1 

0 

0 

c 

1 

0 

0 

1 

0 

1 

0 

1 

0 

0 

1 

D 

Either  of  the  two  remaining  columns  is  acceptable  and  the  row  ZR  is  the 
same,  but  row  2  has  a  zero  value  for  B.   Therefore,  we  choose  either 
(2,2)  or  (2,5).   In  this  case  it  is  immaterial  since  both  are 
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variables.   If  (2,5)  were  chosen  the  result  would  be  to  set  SR(2)  =  kt 
SC(5)  =  kt    R(2,2)  =  -1,  R(2,5)  =  1,  R(5,2)  =  2,  and  R(5,5)  =  0.   The 
final  pivot  is  R(5,2)  meaning  that  its  value  of  2  must  be  used  to 
generate  a  divide  by  2  on  the  fifth  element  of  B.   Since  this  has  not 
been  stored  yet,  we  generate  D/2  ->  B^   as  one  of  the  second  sets  of 

steps.   Note  that  apart  from  the  formation  of  S+T  +  R   ,  no  first  step 
processing  must  be  done  when  the  values  of  L,  R,  S,  and  T  are  changed. 
(This  is  a  fortunate  example  in  this  respect  because  R  is  almost 
triangular  after  permutation.) 

The  back  substitution  step  requires  that  we  handle  the  pivots 
in  the  reverse  order.   Thus,  we  find  the  fifth  pivot  in 
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5 

2 

1 

k 

1 

2 

1 

1 

2 

3 

3 

1 

-1 

0 

0 

X 

1 

U 

2 

0 

-1 

0 

0 

1 

0 

2 

3 

R 

-R 

1 

0 

0 

c 

1 

2 
2 

1 

0 

0 

1 

0 

1 
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0 

1 

0 

0 

0 

B^ 

It  is  the  (5,2)  element.   All  elements  in  the  second  column  must  be 
eliminated  by  generating  the  steps 


1  +  B5  *  Bl 


C  +  R*B   ■>  B 


and  modifying  the  sixth  column  of  R  to  read 
CB1  B5  B3  1  B5]T 
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The  next  pivot  element  is  (2,5)  and  the  process  generates 


B  -  x*B  ->  B 


for  the  one  nonzero  nonpivot  element  in  the  fifth  column. 
Continuing,  we  generate 

B3  -  R*Bl  -  B3 

1  "  Bl  +  BU 

(and  set  B(U)  =  R(U,6)  =  B,  ) 

"to  complete  the  process.   The  reorganization  saved  two  steps  in  SET  1 
and  one  step  in  SET  2.   For  larger  matrices  it  could  be  much  more 
significant.   Note  that  the  solution  obtained  needs  reordering.   B(l) 
to  B(5)  corresponds  to  x(l),  x(5)»  x(3),  x(U),  x(2)  because  not  all 
pivots  were  diagonal  elements. 

There  is  some  choice  of  storage  organization.   For  each 
nonzero  element  it  is  desirable  to  store  its  type  (complexity), 
value  (or  pointer  to  memory  location  containing  a  character  string  o 
address  of  place  in  which  it  will  be  available  at  execution  time), 
and  its  row  and  column  number.   Since  it  is  necessary  to  scan  both 
rows  and  columns ,  a  chain  in  both  row  and  column  directions  is 
desirable.   For  each  row  and  column  we  need  the  nonzero  counts,  we 
need  to  know  if  a  row  has  already  been  handled  and  we  need  to  be 
able  to  trace  backwards  through  the  pivot  elements.   The  latter 
could  be  handled  by  a  chain  through  these  elements,  but  that  would 
make  them  unequal  in  size  from  the  nonpivot  elements,  complicating 
storage  assignment.   Consequently,  it  appears  better  to  construct 
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backward  chains  through  the  row  and  column  vector  elements  which  store 
the  number  of  nonzeros.   (These  are  not  used  after  the  row  and  column 
take  part  in  a  pivot.)   Thus,  the  original  form  of  R  is  a  compacted 
matrix  represented  by  row  and  column  vectors  containing  the  chain 
pointers.   The  result  of  "inverting"  R  will  be  another  compacted  matrix 
whose  elements  may  be  numeric,  symbolic,  or  addresses  of  locations 
containing  the  result.   (In  the  latter  case  code  will  have  been  compiled 
to  produce  the  result.)   The  "lower"  part  of  the  matrix  can  represent 
the  elimination  steps  yet  to  be  performed  (those  that  were  not  either 
produced  as  code  or  performed  numerically),  while  the  diagonal  represents 
the  division  steps  and  the  "upper"  triangular  part  the  back  substitution 
steps  . 

Singularities  will  be  detected  in  R  when  we  find  that  the 
minimum  number  of  nonzero  elements  in  a  remaining  column  is  zero.   That 
column  should  be  ignored  until  the  elemination  process  has  been 
completed.   At  that  time  there  will  be  a  number  of  zero  rows  of  R  left, 
and  a  number  of  columns  not  used  in  the  pivoting  process.   Those  rows 
containing  nonzero  elements  on  the  righthand  side  of  the  reduced  B  can  be 
eliminated.   The  equations  corresponding  to  the  others  can  be  placed  in 
the  nonlinear  set  (3-2)  and  associated  with  variables  in  columns  that  were 
not  used  as  pivots. 

Considerable  theoretical  study  remains  to  be  done  to  determine 
if  this  procedure  will  always  work  if  the  network  is  reasonably  well 
specified  and  to  find  methods  that  will  determine  when  the  network  is 
improperly  specified. 
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5 .   CONCLUSIONS 

A  system  which  is  currently  being  implemented  has  "been 
outlined.   When  complete,  this  system  will  accept  arbitrary  definitions 
of  elements,  and  be  capable  of  doing  transient  analyses  of  correctly 
specified  networks.   Most  incorrect  networks  will  be  detected.   There 
are  many  problems  and  extensions  to  be  analyzed,  and  it  is  planned 
to  continue  the  analysis  throughout  the  implementation  phase.   The 
major  problem  is  that  of  relating  detected  errors  to  the  source  of 
the  error  in  user  notation.   Singularities  in  the  system  of  equations 
could  be  caused  by  any  one  of  many  elements  or  connections,  and  it 
is  necessary  to  tell  the  user  the  most  likely  causes. 

Extensions  that  need  to  be  analyzed  include  the  following. 
Discontinuous  functions  possibly  corresponding  to  the  use  of 
discrete  variables  mixed  and  logical  simulation  problems,  would  allow 
a  useful  class  of  elements  to  be  introduced.   Probabalistic 
functions  would  also  help  in  many  simulations,  although  they  would 
greatly  complicate  the  numerical  algorithms.   Time  delay  elements  are 
important  in  control  systems  and  should  be  allowed,  although  they  also 
complicate  good  numerical  integration  algorithms.   In  more  generality, 
convolution  functions  would  allow  integral  equations  to  be  tackled. 
Finally,  arbitrary  arithmetic  functions  defined  by  subroutines  would 
be  desirable;  these  raise  many  problems  because  the  derivatives  of  the 
functions  are  required  in  the  analysis  presently  used. 
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